Voronoi cell finite element method for heat conduction analysis of composite materials

In this paper, Voronoi cell finite element method (VCFEM) based on assumed flux hybrid formulation has been presented for heat conduction problem of particle reinforced composites material. The heat fluxes satisfying a priori internal thermal balance are directly approximated independently in the matrix and the inclusion respectively. The temperatures on element boundary and matrix-inclusion interface are interpolated by nodal temperature. The thermal balance on the interelement boundary and matrix-inclusion interface is relaxed and introduced into the functional by taking the temperature as Lagrange multiplier. In this way, a functional containing two variables of heat flux and temperature is proposed. Full field heat flux and effective thermal conductivity are obtained. Feasibility and effectiveness of the proposed approach are verified through several numerical examples.

VCFEM 11 based on hybrid flux model and transformation temperature gradient method for materials imbedded with holes were proposed for steady-state heat conduction problem.The variational theory used in the heat conduction process is similar to the transformation strain method in the elastic theory which has proved to be less effective compared with direct constraint method.Special Voronoi elements 12 based on fundamental solutions are built for analyzing clustering on thermal effect in fiber reinforced cement composites.Fundamental solutions need to be constructed for special example, it does not have universal applicability for arbitrary shapes and fibers shape in that paper are all circler.
Considering advantages of VCFEM in calculating particle reinforced composites, a new VCFEM based on hybrid flux model and direct constraint method is proposed for heat conduction problem.In this model, the heat flux in the domain of each element and the temperature on the element boundary are assumed independently.The heat flux in the matrix and the inclusion satisfying thermal balance are assumed independently.Thermal balance on interelement boundary and matrix-inclusion interface inside the element are relaxed by temperature as Lagrange multiplier introduced into the functional.Thus, the functional containing two variables of heat flux and temperature is built ("Functional derivation").Heat flux approximations satisfying equilibrium relations in matrix and inclusion are built up from the heat flux function respectively ("Equilibrated flux fields in VCFEM").Temperatures on element boundary and matrix-inclusion interface are interpolated by nodal temperature.The heat flux parameters inside the element can be linked with the displacements of the internal node and the external node of the element, and can be condensed together with the displacements of the internal node, thus obtaining the same standard solution form as the traditional FEM ("Method of solution").Numerical examples are shown in "Numerical examples" and conclusions are drawn in "Conclusions".

Functional derivation
Microstructure containing inclusions is tessellated by Voronoi mesh, with each element containing an inclusion as shown in Fig. 1b.The matrix region of each element is Ω m , the inclusion region is Ω c .In the Fig. 1a it is seen that the element boundary Γ consists of inter-element boundary Γ e , prescribed temperature boundary Γ θ and prescribed heat flux boundary Γ q with outward normal n e (Γ = Γ e ∪ Γ θ ∪ Γ q ).Γ c is matrix-inclusion interface with an outward normal n c = n c x ,n c y .The heat flux vector q m with cartesian component ( q m x ,q m y ) in the Ω m and the heat flux vector q c with cartesian component ( q c x ,q c y ) in the Ω c are assumed based on satisfying the self-equilibrating heat flux as shown in Eq. (1) (excluding external heat transfer).This will be explained and justified in "Equilibrated flux fields in VCFEM".And Fourier's law for the constitutive relation of heat flux is: where k m and k c are the thermal conductivity matrices of the matrix and inclusion, respectively.θ m and θ c are the corresponding temperatures.
(1)  where θ e , θ and q n are the displacement of the element boundary in Eq. ( 9), the prescribed temperature on Γ θ and the prescribed heat flux on Γ q , respectively.The outward normal heat flux q m n of the element and the heat flux q c n in the outward normal direction at the matrix-inclusion interface can be expressed as As shown in Fig. 1c, the outer normals n e , n m and n c at the inter-element boundary and matrix-inclusion boundary have the following transformation relations: where the superscripted plus and minus signs are just to visually show the orientation of the normals of neighboring regions on the common boundary.
Further, the following conditions of heat flux continuity are also satisfied on the matrix-inclution boundary Γ c and the inter-element boundary Γ e .
where q + and q -denote the heat fluxes of neighboring cells on the inter-cell boundary, respectively.
The energy functional e may be expressed for each element in terms of heat flux and boundary/interface temperature fields as Equation (8) can also be expressed as: where Γ = Γ e ∪ Γ θ ∪ Γ q , and the element boundary temperature θ e is equal to θ when located at the prescribed temperature boundaries Γ θ .
The total energy for the entire domain is obtained by adding each element functional e as: S m and S c are the inverse of isotropic thermal conductivity matrix that are expressed as: Compatible temperature fields are continuous on the inter-element boundary and the matrix-inclusion interface.The temperature θ e at the element boundary and the temperature θ c at the matrix-inclusion interface are obtained by interpolating the nodal temperatures of the following form: where L is the interpolating function, e and c are temperature at node of element boundary and matrix- inclusion boundary respectively.
(3) θ e =θonŴ θ (4) n e q m = q n on Ŵ q (5) n =n e q m on Ŵ e q m n = n m q m on Ŵ c q c n =n c q c on Ŵ c (6) n e =n e+ on Ŵ e n c =n c + ,n m =n c -on Ŵ c (7) n c− q m +n c + q c = 0 on Ŵ c n e + q + +n e− q − = 0 on Ŵ e (8) n e+ q + + n e− q − θ e dŴ + Ŵ θ n e q m θ dŴ + Ŵ q (n e q m − qn )θ e dŴ � e = 1 2 �m q mT S m q m d�+ 1 2 �c q cT S c q c d� + Ŵ n e q m θ e dŴ − Ŵq q n θ e dŴ + Ŵc n c− q m + n c+ q c θ c dŴ www.nature.com/scientificreports/

Equilibrated flux fields in VCFEM
In VCFEM formulation, the equilibrium conditions Eq. (1) in the matrix and inclusion respectively are satisfied a priori in a strong sense.In the absence of heat source, heat flux fields satisfying equilibrium relations can be built up from the heat flux function .In the local coordinate system, the resulting heat flux are where ε and η are localized coordinate as in Fig. 2. Since bringing Eq. ( 13) into Eq.(1) yields a constant relationship, we assume that the heat flux field assumed in the form of Eq. ( 13) is presupposed to satisfy the equilibrium conditions.
Different functional forms of the heat flux functions are chosen for the matrix and inclusion phases in VCFEM implementation.Independent construction of heat flux functions m and c in matrix and inclusion allow for flux discontinuities across the matrix-inclusion interface.
The proper selection of heat flux function has a significant influence on the convergence and efficiency of VCFEM.The construction of heat flux functions should consider the influence of shape location and distribution of inclusion.The influence on the shape of matrix heat flux functions far from the interface should disappear.Heat flux function for matrix is constructed of pure polynomial forms m poly and reciprocal functions 13 m rec , and heat flux function for inclusion is constructed of pure polynomial forms c poly .The heat flux functions for the matrix and inclusion are constructed as: The local coordinate system takes the center of the inclusion as the origin, and the horizontal and longitudinal axes correspond to the long half axis and the short half axis of the inclusion respectively.In the above function, the pure polynomial part is constructed from the Pascal triangle in the local coordinates (ε , η) shown in Fig. 2, which is given by the following formula: The term of the reciprocal function is expressed as: The function f has the following property when (ε, η) → 0 , 1/f → 0 , and can be used as a mapping function to construct a shape-based reciprocal heat flux function.Furthermore, at the matrix-inclusion interface, f = 1, and the coefficient �β is equal to the polynomial heat flux coefficient, ensuring that heat flux is continuous at the interface.In the region away from the inclusions, the effect of the interaction heat flux function can be overlooked.Thus, near the matrix-inclusion interface, the heat flux is a combination of polynomial and interaction components, and the closer to the interface, the stronger the effect, ensuring that the approach can better compute the heat flux concentration at the interface.According to this characteristic, it can be taken as an elliptic function.
The associative Eqs. ( 13)-( 16), the heat flow density can be expressed as: www.nature.com/scientificreports/Therefore, each indeterminate field is expressed as product of the function P of localized coordinate alone and the heat flux coefficients β alone: Take Eqs. ( 12) and (19) into Eq.( 9): where the H-matrix and G-matrix are defined as follows: .
Where Γ = Γ e ∪ Γ θ ∪ Γ q , and the element boundary temperature θ e is equal to θ on the prescribed temperature boundaries Γ θ .'m' , and 'c' are superscripts used to distinguish the matrix and inclusion parts, respectively.

Method of solution
Setting the primary variations of e in Eq. (20) with respect to the heat flux coefficients β m and β c , respectively, to zero: yields Heat flux coefficients β from Eq. ( 23) can be written as Setting the primary variations of the total energy functional in Eq. ( 10) with respect to e and c to zero: gives Take Eq. (24) into Eq.( 26) The second part relates only within each element, since the internal node temperature of each element is only related to the external node temperature of itself.The relationship between them is satisfied as follows: Substituting Eq. (31) into Eq.( 29), the solved relation can be expressed as: where, K e =K 11 − K 12 K −1 22 K 12 is the element stiffness matrix.The computational structure of the program mainly consists of an input module, an initial computation module, a solution module and an output module.Its general calculation flow is shown in Fig. 3.

Numerical examples A convergence study for heat flux
A convergence study of one element A study is conducted to examine the sensitivity and convergence of the VCFEM.A square plate with side length of 100 mm containing an inclusion of radius 15 mm is studied.Equations ( 15) and ( 16) with various numbers of terms will be utilized for the computation to strike a balance between calculation efficiency and correctness.When the number of terms for Eqs. ( 15) and ( 16) is taken as 20 and 14, respectively, it is discovered that the VCFEM has greater computing efficiency with guaranteed correctness.Thus, the heat flux field in the inclusion is generated with a complete 6th order heat flux function Eq. (15) corresponding to 20 β terms in the flux poly- nomial.The matrix flux field has an additional 14 reciprocal terms due to the 5 th reciprocal heat flux function in Eq. ( 16).This result a total of 34 β in matrix.The thermal conductivity k m is 180 W/(mK), and the thermal conductivity k c is 330 W/ (mK).Temperature of 10 K was applied to all nodes on the downside and temperature of 0 K was applied to all nodes on the upside.
The results of the VCFEM calculations were compared with those of the finite element software MARC.For MARC, two groups of meshes are established, one is a coarse mesh with 3418 triangular elements, the other is a fine mesh with 7141 triangular elements.All the boundary conditions are the same of Voronoi mesh.A comparison of heat flux in y direction at mid-section) is shown in Fig. 4. It can be seen from the figure that the results calculated using a large number of 7141 finite elements are agreed well to those calculated by using VCFEM.However, the results calculated using fewer finite element of 3418 elements deviate from the accurate values.This proves high computational efficiency and accuracy of the proposed VCFEM.
Furthermore, in order to verify the applicability of VCFEM to generally shaped materials, a non-rectangular shaped model was computed with VCFEM (Fig. 5d-f) and the results were compared with those of the finite element software MARC(Fig.5a-c).Where 34β is still used for the matrix and 20 polynomial terms are used for inclusions.The cloud diagram of the heat flux is shown in Figs. 5, and 6 illustrates the comparison results on the sampling path in Fig. 5b.The results show that the VCFEM calculations are generally satisfactory for models with pointed/cornered profiles, but the accuracy at the pointed corners needs to be improved as shown at the right end of Fig. 6b.Sensitivity to geometric distortion This section's primary goal is to examine, using the model with four Voronoi cells in Fig. 7, how sensitive Voronoi cells are to geometric distortions.The model's upper and lower boundaries have temperatures of 0 K and 10 K, respectively.The left and right boundaries are adiabatic, and the matrix and inclusions' thermal conductivities are 180 W/mK and 330 W/mK, respectively.Figure 7a shows the most ideal case, and the cells of Fig. 7b, c have been partitioned into the grid using linear and parabolic forms, respectively 14 .In the calculations, it is assumed that the values are positive when the resulting rightward and upward δ.The geometric distortion is determined by the parameter 2δ/l, which fluctuates in the range −0.3:0.3 to ignore as much as possible the effect due to the close distance between the inclusion and  the cell boundary.Figure 7 shows how the heat flux varies with 2δ/l for the matrix and inclusions at A and B, respectively.Calculations in Fig. 8 demonstrate that when the mesh produces geometrical distortions, changes on the order of 10 -2 are produced in the heat fluxes of the matrix and inclusions.As a result, the Voronoi cell's effect on geometrical distortion when calculating heat flux can be judged to be within an acceptable error range.This suggests that the effects of geometric distortion in conventional finite element methods have been somewhat mitigated by VCFEM.

A convergence study for RVE containing multiple inclusions
In this example, RVE (2.5 mm × 2.5 mm) containing 16 elliptical inclusions with size and location randomly distributed is considered.The volume fraction of inclusions is 5%.Figure 9 shows mesh by MARC and VCFEM models.The MARC mesh consists of 14,860 Tri elements while the VCFEM mesh has only 16 elements corresponding to the number of inclusions.The thermal conductivity of the matrix is 180W/mK and the thermal conductivity of inclusions is 330W/mK.A temperature of 10 K is applied to the bottom edge and a temperature of 0 K is applied to the top edge of the model.Figure 10 shows the true microscopic heat flux at the midsection(y/L = 0.5).Full field heat flux distribution is shown in Fig. 11.

Effect of microstructure on response
In this section the effect of volume fraction, shape, thermal conductivity and number of inclusions, on the macroscopic as well as true microstructural response of RVEs is studied.Representative volume elements (RVE) characterize a point by microstructure.The effective heat conductive of RVE correspond to average heat flux divided by the temperature gradient.The macroscopic average heat flux corresponds to the volume average of the microscopic heat flux.When apply fixed temperature 2 on the upside and 1 on the downside, the vertical average heat flux of the RVE is calculated as  www.nature.com/scientificreports/where, L is the side length of a square RVM.The temperature gradient in the Y direction was calculated from: The Voronoi cell finite element model uses 34 β parameters in the heat flux interpolation for matrix and 20 β parameters in the heat flux interpolation for inclusion for all case parameters in the stress interpolation for all cases.

Effect of volume fraction of inclusion
In this study, the thermal conductivity k m of matrix is 180 W/mK, and the thermal conductivity k c of reinforcement is 330 m W/mK.Various volume fractions of circular inclusions are considered.The macroscopic response effective thermal conductivity (k e ) increases with the increase of inclusion volume fraction, as shown in the Fig. 12.This is expected since larger volume of the inclusions will cause the overall response to move in the direction of its individual response.At the same time, the analytical value calculated by flexible model 12 as Eq. ( 35) is also drawn in the Fig. 12.When the volume fraction of inclusion is 5%, the relative error between the value calculated by VCFEM and the analytical value is 0.1%, when the volume fraction of inclusion is 40%, the relative error between the value calculated by VCFEM and the analytical value is 1.4%.Figure 13 shows the y-direction heat flux cloud map corresponding to the volume fraction of the inclusions mentioned above.

Effect of shape, orientation and thermal conductivity of inclusion
In this study, the thermal conductivity k m of matrix is 180 W/(m K), and the thermal conductivity k c of reinforcement varies from 30 to 480 W/(m K) and increases gradually at intervals of 50 W/(mK).The effect of shape and orientation together with of thermal conductivity of inclusions on the behavior of microstructure and macrostructure is investigated respectively.Sixteen elliptical inclusions of various shapes, orientation and thermal conductivity but with a constant volume fraction of inclusion of 5%, are distributed randomly in the RVE. Figure 9b shows the VCFEM mesh for the fifth group, and Fig. 15 shows the VCFEM meshes for the other four different distributions.They are (1) major axis in the horizontal direction (2) major axis in the vertical direction (3) major axis rotates 45° from the horizontal direction (4) major axis rotates 145° from the horizontal direction (5) arbitrary shapes, sizes and orientation.The ratio of major to the minor axis, for ellipses in the last four cases is 2.3.For each of the five groups, the effective thermal conductivity changing with different thermal conductivity of particles were calculated and shown in Fig. 2. It can be seen from the Fig. 2 that when the thermal conductivity of the particles increases from 30W/mK to 270 W/mK, k e of the composite material shows an increasing trend.however, the effects of particle shape on the effective thermal conductivity are different with different particle thermal conductivity.
(1) When 30W/mK ≤ k c < 130 W/mK, the shape and direction of the particles have great influence on the macroscopic k e as shown in Fig. 14.At this interval, the k m is higher, the K c is lower, and the matrix is the main carrier for conducting heat.When the orientation of the particles is different and the distance between the edges of the particle's changes, the width and path of the heat conduction channel of the matrix will change.When the particles are distributed in the group 1, the distance between the edges of the particles is small, thus the width of the heat conduction channel of the matrix becomes narrow, and the matrix is hindered more on the heat conduction path.Therefore, k e of the composite material is reduced, and k e of microstructures of the group 2 is 162.69W/mK which is the smallest among the five groups, and k e of the group 3 is 170.16W/mK which is the largest.The difference between the maximum value and the minimum value is 7.47W/mK.k e of other microstructures is between therm.Taking k m = 180 W/mK, k c = 30 W/mK for example, Fig. 15 is a contour of microscopic heat flux distribution of composite materials for group 1 and group 2 simulated by VCFEM.From Fig. 15 we can see that when the heat flux meets the particles   It can therefore be conducted that microscopic as well as macroscopic response is highly sensitive to the distribution of shapes of inclusions when k c is deviated from k m , while those response is insensitive to the distribution of shapes of inclusions when k c is equivalent with k m.

Effect of inclusion numbers
In this study, the effect of number of inclusions on the behavior of microstructure and macrostructure is investigated.Inclusions of various shapes, sizes and orientation are distributed randomly in the RVE.Four groups of different inclusion numbers with a constant volume fraction of V c = 5%, are studied.These four groups contain50, 100, 500 and 16 inclusions respectively.Matrix k m = 180 W/mK and inclusion k c = 330 W/mK are used for this study.Contour plot of microstructures with different number of 50, 100 and 500 inclusions are shown in Fig. 16, while 16 inclusions is shown in Fig. 16b.From these figures, it can be seen that the heat flux distribution of microstructures containing 16, 50 and 100 inclusions is very similar, the difference between the maximum value is very

Figure 1 .
Figure 1.(a) The typical Voronoi element.(b) Voronoi mesh for materials with various inclusions.(c) The schematic diagram of each outer normal of 'element A' as the research object.

Figure 4 .
Figure 4. Flux distribution along x/L = 0.5.(a) Comparison of the results of different combinations of Poly and Rec terms with MARC calculations; (b) the results of VCFEM and MARC with different number of cells are calculated.

Figure 5 .
Figure 5. Mesh delineation and computational results for non-rectangular models.(a,d) MARC's meshing scheme with 4754 elements and a VCFEM cell; (b,e) horizontal heat flux q x calculated by MARC and VCFEM; (c,f) vertical heat flux q y calculated by MARC and VCFEM.

Figure 6 .
Figure 6.Heat flux on the sampling path in Fig. 5. (a) Horizontal heat flux q x comparison results, (b) vertical heat flux q y comparison results.

Figure 7 .
Figure 7. Meshing scheme for Voronoi elements: (a) normal Voronoi cells, (b) Voronoi cells with small linear geometric distortions and (c) Voronoi cells with small parabolic geometric distortions.

Figure 8 .
Figure 8.(a) Sensitivity of the heat flux density of the matrix at A to geometric distortion.(b) Sensitivity of the heat flux density of inclusions at B to geometric distortion.

Figure 12 .
Figure 12. k e computed by VCFEM and analytical method.

Figure 13 .Figure 14 .
Figure 13.Heat flux fields in y direction with volume fracture of inclusion varies from 5 to 40% at 5% intervals.